数学原理简介¶
多元回归中的普通最小二乘法是拟合参数向量 $\beta$, 使得 $\|X \beta-Y\|_{2}^{2}$ 达到最 小值。岭回归是选择合适的参数 $k \geq 0$, 拟合参数向量 $\boldsymbol{\beta}$, 使得 $\|X \boldsymbol{\beta}-\boldsymbol{Y}\|_{2}^{2}+\boldsymbol{k}\|\boldsymbol{\beta}\|_{2}^{2}$ 达到最小值, 解决了 $X^{T} X$ 不可逆的问题。LASSO 回归, 是 选择合适的参数 $k \geq 0$, 拟合参数向量 $\beta$, 使得 $$ J(\beta)=\|X \beta-Y\|_{2}^{2}+k\|\beta\|_{1} $$ 达到最小值，其中 $\boldsymbol{k}\|\boldsymbol{\beta}\|_{1}$ 为目标函数的惩罚项, $\boldsymbol{k}$ 为惩罚系数。
由于式（12.17）中的惩罚项是关于回归系数 的绝对值之和，因此惩罚项在零点处是不可导的，那么应用在岭回归上的最小二乘法在此失效，不仅如此，梯度下降法、牛顿法与拟牛顿法都无法计算出LASSO回归的拟合系数。坐标下降法可以求得LASSO回归系数，坐标下降法与梯度下降法类似，都属于迭代算法，所不同的是坐标轴下降法是沿着坐标轴下降，而梯度下降则是沿着梯度的负方向下降，具体的数学原理这里就不介绍了。
由于拟合LASSO回归模型参数时，使用的损失函数（机器学习中的用语）式（12.17）中包含惩罚系数 ，因此在计算模型回归系数之前，仍然需要得到最理想的$k$值。与岭回归模型类似，$k$值的确定可以通过定性的可视化方法。
案例¶
$$\begin{array}{ccccc|ccccc} \hline \text { 年份 } & x_{1} & x_{2} & x_{3} & y & \text { 年份 } & x_{1} & x_{2} & x_{3} & y \\ \hline 1949 & 149.3 & 4.2 & 108.1 & 15.9 & 1955 & 202.1 & 2.1 & 146.0 & 22.7 \\ 1950 & 171.5 & 4.1 & 114.8 & 16.4 & 1956 & 212.4 & 5.6 & 154.1 & 26.5 \\ 1951 & 175.5 & 3.1 & 123.2 & 19.0 & 1957 & 226.1 & 5.0 & 162.3 & 28.1 \\ 1952 & 180.8 & 3.1 & 126.9 & 19.1 & 1958 & 231.9 & 5.1 & 164.3 & 27.6 \\ 1953 & 190.7 & 1.1 & 132.1 & 18.8 & 1959 & 239.0 & 0.7 & 167.6 & 26.3 \\ 1954 & 202.1 & 2.2 & 137.7 & 20.4 & & & & & \\ \hline \end{array}$$下表 是 Malinvand 于 1966 年提出的研究法国经济问题的 一组数据。所考虑的因变量为进口总额 $y$, 三个解释变量分别为: 国内总产 值 $x_{1}$, 储存量 $x_{2}$, 总消费量 $x_{3}$ (单位均为 10 亿法郎)。
求例 $12.3$ 的 LASSO 回归方程。
解 画出的 $k$ 与LASSO回归系数的关系图如图12.2所示, 从图中可以 看出选 $\boldsymbol{k}=\mathbf{0 . 2 1}$ 较好。对应的标准化LASSO回归方程为 $$ \hat{y}^{*}=0.0136 x_{2}^{*}+0.7614 x_{3}^{*}, $$ 将标准化回归方程还原后得 $$ \hat{y}=-1.6602+0.0374 x_{2}+0.1677 x_{3} \text {, } $$ 模型的拟合优度 $R^{2}=0.9061$ 。从计算结果可以看出, LASSO 回归可以非常 方便地实现自变量的笑选。
代码
import numpy as np; import matplotlib.pyplot as plt
from sklearn.linear_model import Lasso, LassoCV
from scipy.stats import zscore
plt.rc('font',size=16)
plt.rc('text', usetex=True) #没装LaTeX宏包把该句注释
a=np.loadtxt("data/economic.txt")
n=a.shape[1]-1 #自变量的总个数
aa=zscore(a) #数据标准化
x=aa[:,:n]; y=aa[:,n] #提出自变量和因变量观测值矩阵
b=[] #用于存储回归系数的空列表
kk=np.logspace(-4,0,100) #循环迭代的不同k值
for k in kk:
md=Lasso(alpha=k).fit(x,y)
b.append(md.coef_)
st=['s-r','*-k','p-b'] #下面画图的控制字符串
for i in range(3): plt.plot(kk,np.array(b)[:,i],st[i]);
plt.legend(['$x_1$','$x_2$','$x_3$'],fontsize=15); plt.show()
mdcv=LassoCV(alphas=np.logspace(-4,0,100)).fit(x,y);
print("最优alpha=",mdcv.alpha_)
#md0=Lasso(mdcv.alpha_).fit(x,y) #构建并拟合模型
md0=Lasso(0.21).fit(x,y) #构建并拟合模型
cs0=md0.coef_ #提出标准化数据的回归系数b1,b2,b3
print("标准化数据的所有回归系数为：",cs0)
mu=np.mean(a,axis=0);
s=np.std(a,axis=0,ddof=1) #计算所有指标的均值和标准差
params=[mu[-1]-s[-1]*sum(cs0*mu[:-1]/s[:-1]),s[-1]*cs0/s[:-1]]
print("原数据的回归系数为：",params)
print("拟合优度：",md0.score(x,y))
案例2¶
\begin{array}{cccccc} \hline \text { 年 } & x_{1} & x_{2} & x_{3} & x_{4} & y \\ \hline 1994 & 3496.2 & 3.43 & 136.69 & 111.78 & 205.42 \\ 1995 & 4283 & 3.52 & 145.27 & 115.7 & 249.96 \\ 1996 & 4838.9 & 3.73 & 147.52 & 118.58 & 289.67 \\ 1997 & 5160.3 & 3.94 & 158.25 & 122.64 & 358.36 \\ 1998 & 5425.1 & 4.16 & 163 & 127.85 & 423.65 \\ 1999 & 5854 & 4.37 & 183.2 & 135.17 & 533.88 \\ 2000 & 6280 & 4.59 & 207 & 140.27 & 625.33 \\ 2001 & 6859.6 & 4.81 & 234.17 & 169.8 & 770.78 \\ 2002 & 7702.8 & 5.02 & 325.1 & 176.52 & 968.98 \\ \hline \end{array}在建立中国私人轿车拥有量模型时, 主要考虑以下因素: (1) 城镇居民家庭人均可支配收入 $x_{1}$ (元)，(2) 全国城镇人口 $x_{2}$ (亿人), (3) 全国汽车产量 $x_{3}$ (万辆), (4) 全国公路长度 $x_{4}$ (万公里) 。具体数据见表 $12.3$, 其中 $y$ 表示中国私人轿车拥有量 (万䢀)。试建立 $y$ 的经验公式。
对于上述问题, 我们可以直接用普通的最小二乘法建立 $y$ 关于四个解释 变量 $x_{1}, x_{2}, x_{3}$ 和 $x_{4}$ 的回归方程为 $$ \hat{y}=-1028.4134-0.0159 x_{1}+245.6120 x_{2}+1.6316 x_{3}+2.0294 x_{4} \text {, } $$ 模型的检验见下面程序运行结果, 我们这里就不具体给出了。 在 $\alpha=0.05$ 的水平下, 以上的回归方程是显著的。但变量 $x_{1}$ 对 $y$ 是不显 著的, 且回归方程中 $x_{1}$ 前面的系数为负值也不合理。
我们选择 $k=0.05$, 建立的 LASSO 回归方程为 $$ \hat{y}=-908.2059+203.0938 x_{2}+1.4562 x_{3}+2.0469 x_{4} \text {. } $$
import numpy as np; import matplotlib.pyplot as plt
import statsmodels.api as sm
from sklearn.linear_model import Lasso
from scipy.stats import zscore
#plt.rc('text', usetex=True) #没装LaTeX宏包把该句注释
a=np.loadtxt("data/car.txt") #加载表中的9行5列数据
n=a.shape[1]-1 #自变量的总个数
x=a[:,:n] #提出自变量观测值矩阵
X = sm.add_constant(x)
md=sm.OLS(a[:,n],X).fit() #构建并拟合模型
print(md.summary()) #输出模型的所有结果
aa=zscore(a) #数据标准化
x=aa[:,:n]; y=aa[:,n] #提出自变量和因变量观测值矩阵
b=[] #用于存储回归系数的空列表
kk=np.logspace(-4,0,100) #循环迭代的不同k值
for k in kk:
md=Lasso(alpha=k).fit(x,y)
b.append(md.coef_)
st=['s-r','*-k','p-b','^-y'] #下面画图的控制字符串
for i in range(n): plt.plot(kk,np.array(b)[:,i],st[i]);
plt.legend(['$x_1$','$x_2$','$x_3$','$x_4$'],fontsize=15); plt.show()
md0=Lasso(0.05).fit(x,y) #构建并拟合模型
cs0=md0.coef_ #提出标准化数据的回归系数b1,b2,b3,b4
print("标准化数据的所有回归系数为：",cs0)
mu=a.mean(axis=0); s=a.std(axis=0,ddof=1) #计算所有指标的均值和标准差
params=[mu[-1]-s[-1]*sum(cs0*mu[:-1]/s[:-1]),s[-1]*cs0/s[:-1]]
print("原数据的回归系数为：",params)
print("拟合优度：",md0.score(x,y))
